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INTRODUCTION 

Programs  for  non-linear  parameter  estimation  exist  for  a  variety  of  eomputer  languages, 
many  of  whieh  are  not  eompatible  with  eurrent  PCs,  and  most  programs  are  not  user  friendly. 

To  make  a  program  that  ean  be  run  on  a  PC,  by  users  unfamiliar  with  eomputer  programming,  a 
modified  non-linear  parameter  estimation  program  was  written  for  Matlab.  The  program  allows 
parameter  estimation  to  be  made  by  the  method  of  least  squares,  and  a  modifieation  allows 
estimation  based  on  the  method  of  maximum  likelihood  (1). 

The  non-linear  parameter  estimation  method  is  based  on  the  approaeh  by  Marquardt  (5), 
with  a  modifieation  allowing  maximum  likelihood  estimation  (1).  Briefly,  it  ean  be  shown  that  if 
a  parameter  Lambda  is  ehosen  to  be  large  enough,  the  parameters  (P)  will  always  eonverge  at  the 
value  giving  the  best  fit  by  the  least  squares  eriterion  (5).  The  smaller  the  value  of  Lambda,  the 
faster  the  program  will  reaeh  eonvergenee  (5).  For  an  initial  Lambda  and  a  set  of  starting  P 
values,  the  program  will  ealeulate  a  new  set  of  P  values.  Next,  the  program  eompares  the  sum  of 
squared  errors  (SSL)  or  the  log-likelihood  (LL)  for  the  eurrent  set  of  parameters  and  eompares  it 
with  the  SSL  or  LL  for  the  old  P  values.  If  the  new  set  of  P  values  reduees  the  deviation 
(reduees  the  SSL  or  inereases  the  LL),  a  new  set  of  P  values  are  eomputed  by  first  redueing  the 
Lambda  by  a  faetor  of  10.  Conversely,  if  the  deviation  is  inereased,  the  new  set  of  P  values  are 
eomputed  by  first  inereasing  the  Lambda  by  a  faetor  of  10.  This  iterative  proeedure  eontinues 
until  the  program  has  found  a  new  set  of  P  values  that  does  not  ehange  the  devianee  by  more 
than  a  set  value  from  the  old  P  values,  the  eonvergenee  eriterion  (eonv).  When  this  is  aehieved, 
the  program  eomputes  the  result  for  the  eonverged  set  of  parameters  using  propagation  of  error 
formulas  (4).  For  mathematieal  justifieation  for  the  routines,  the  interested  reader  is  referred  to 
the  referenees  (1,  3,  5). 

THE  PROGRAM 

The  parameter  estimation  routine,  ealled  “Marquardt”,  was  written  in  the  Matlab 
language  (Matlab  Student  Edition  Version  5).  The  program  is  detailed  in  Appendix  A,  and  the 
routines  required  to  run  the  program  are:  marquardt,  modl-mod5,  E,  get  data,  get  funetion, 
new  array,  binary  and  a  funetion  in  the  form  AB  Homerl .  The  program  uses  funetions  instead 
of  “goto”  statements,  the  advantage  being  a  more  modular  program  that  is  easier  to  read  and 
modify  (6). 
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RUNNING  THE  PROGRAM 

Save  the  data  in  a  text  file,  where  dependent  (y)  and  independent  (x)  variables  are  ordered 
in  eolumns,  starting  with  the  dependent  variable  as  shown  in  Table  1  for  the  data  set  in  Example 
1. 

Before  running  the  program,  the  number  of  independent  variables  (N2),  eonvergenee 
eriterion  (CONV),  number  of  iterations  (Tl),  and  Lambda  (L)  need  to  be  speeified  during 
initialization  in  Marquardt.  Onee  this  is  done,  the  program  is  run  by  typing  “marquardt”  at  the 
“edu>”  prompt  in  Matlab. 

A  diary  file  is  ereated,  named  with  the  date  and  time  and  with  extension  “dry”.  The  diary 
file  saves  the  information  printed  to  the  sereen  in  a  file  that  may  be  opened  with  any  program 
that  is  able  to  open  text  files. 

Next,  the  user  is  asked  to  enter  the  data  file  and  the  funetion  file  name.  The  funetion  file 
eontains  the  number  of  parameters  to  be  used.  Currently,  the  funetion  files  are  named  with  the 
prefix  “AB_”,  but  this  ean  be  ehanged  in  the  funetion  “Get_Funotion”.  After  the  name  of  the 
funetion  has  been  speeified,  the  user  is  asked  to  enter  starting  values  for  the  parameter 
estimation.  Onee  the  starting  values  have  been  entered,  the  parameter  seareh  begins,  with  the 
log-likelihood  or  standard  error  for  eaeh  iteration  printed  to  the  sereen  and  saved  to  the  diary  file. 
At  eonvergenee,  the  final  result,  ineluding  the  parameter  estimates,  the  standard  error  of  the 
parameter  estimates,  the  eoeffieients  of  variation,  and  the  varianee-eovarianee  matrix,  is  printed 
to  the  sereen  and  saved  in  the  diary  file. 

Changes  in  the  program  ean  be  made  for  the  number  of  iterations  (Tl),  the  eonvergenee 
eriterion  (CONV),  the  lambda  (L),  and  the  number  of  independent  variables  (N2).  This  is  done 
by  ehanging  the  initilization  of  any  of  these  in  the  “marquardt”  routine. 

The  program  uses  the  funetion  “binary”  to  determine  if  the  dependent  variable  is  binary 
(1  or  0)  or  eontinuous.  Therefore,  no  ehange  is  neeessary  when  dealing  with  different  data  sets. 
The  proeedure  used  is  seen  in  the  output  of  the  result,  where  eontinuous  data  output  a  SSE,  and 
binary  data  a  EL. 

Addition  of  “mod8”  eomputes  the  95%  eonfidenee  regions  of  the  parameter  estimate. 
Currently,  “mod8”  is  only  able  to  do  this  for  estimation  of  one  independent  variable.  After  the 
estimation  is  eompleted,  “mod8”  plots  the  result.  The  approximation  for  the  standard  error  of  the 
parameters  is  eomputed  using  the  formulas  for  propagation  of  error  (1,4). 
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EXAMPLES 

Example  1,  Continuous  variables 

For  the  estimation  of  eontinuous  variables,  the  data  set  presented  in  Appendix  C  of 
Bailey  and  Homer  (1977)  will  be  used  to  present  the  non-linear  least  square  estimation 
proeedure.  These  data  are  presented  in  Table  1,  and  Appendix  B  shows  the  fit  using  the  non¬ 
linear  equation 

f(x)  =  pi.x/''2*P2^  [1] 

to  the  eontinuous  data. 

Example  2,  Binary  variables 

Many  times  in  medieal  or  biologieal  researeh,  one  is  faeed  with  data  that  are  binary,  i.e. 
response  or  no  response,  death  or  no  death  ete.  For  these  problems,  one  ean  use  a  probabilistie 
formulation.  The  unknown  parameter  beeomes  the  probability  (P)  of  a  speeifie  outeome  or 
response.  Aeeordingly,  the  probability  of  no-response  is  then  1-P. 

In  Appendix  C,  the  binary  data  presented  in  the  Table  2  is  fitted  to  the  dose  response 
funetion. 


II 

• 

[2] 

P(x)=XiP^(piP"  +  pi.Xi)-' 

[3] 

using  the  maximum  likelihood  teehnique.  In  these  eases,  the  outeome  (Y)  is  defined  as 
either  response  (I)  or  no-response  (0).  The  likelihood  of  an  event  for  the  n**'  observation  is: 

L(n)  =  P  .  (1-P)^'-Y("» 

That  is,  the  response  (1)  oeeurs  with  a  probability  P,  and  the  no-response  with  a  probability  1-P. 


The  likelihood  for  the  n  independent  observations  is  the  produet  of  their  outeomes: 

LL=iL(i) 

i=i 

The  estimation  routine  adjusts  the  parameters,  and  eonsequently  P,  to  maximize  the  LL  whieh  is 
defined  as  the  best  fit  to  the  data. 
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Table  1,  Data  format  for  a  file.  The  first  row  shown  should  be  omitted, 

data  should  be  presented  in  the  form  below. 

Reproduced  from  Bailey  and  Homer  (1) 

i,e,  only  the  actual 

Y 

XI 

X2 

1 

0.0001 

0 

0 

0.0001 

2 

4 

1 

1 

2 

1 

2 

8 

2 

1 

2 

2 

0 

Equation  1  is  used  to  fit  the  data 

Table  2,  Binary  data  for  use  in  fitting  with  Equation  2  and  Equation  3, 

Y 

X 

0 

0.6 

0 

0.6 

0 

0.6 

1 

1.04 

1 

1.04 

1 

1.04 

1 

1.44 

1 

1.44 

1 

1.44 

0 

2 

1 

2 

1 

2 

1 

2.75 

1 

2.75 
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APPENDIX  A 

Each  routine  is  separated  by  “%%%%%%%”,  and  its  name  is  written  within  the  string  of  %- 

signs.  In  the  program,  comments  begin  with  a  ”%”-sign,  which  stops  Matlab  from  reading  the 

remainder  of  that  line.  Some  comments  have  been  left  in  for  elarification.  To  use  “modS”,  erase 

the  %-sign  in  front  of  these  in  “marquardf save  the  file  and  run  the  program  again. 

%%%%%%%%%%%%%%BEGINMARQUARDT%%%%%%%%%%%%%% 
function  marquardt ( ) 
clear  all; 

format  long; 

timearray  =  clock; 

month  =  int2str (timearray (2 )) ; 

if  length (month) <2 

month  =  ['O'  month]; 

end 

day  =  int2str (timearray (3) ) ; 
if  length (day) <2 
day  =  ['O'  day] ; 

end 

hour  =  int2str (timearray (4) ) ; 
if  length (hour ) <2 

hour  =  ['O'  hour] ; 

end 

minute  =  int2str (timearray (5) ) ; 
if  length (minute) <2 

minute  =  ['O'  minute]; 

end 

diaryfile  =  [month  day  '-'  hour  minute  '.dry']; 
diary (diaryfile) ; 

disp(['Open  diary  file:  '  diaryfile]); 

%  read  data  from  text  file 

%  input  file  has  records  in  the  following  format  {Y,X(1),  %(2),...X(k) 

%only  one  independent  variable  is  allowed 

%We've  defined  arrays  and  vectors  that  will  change  and  which  needs 

%to  be  accesed  through  the  whole  program  as  global  variables,  while 

%the  other  variables  that  must  not  be  changed  are  passed  to  each  %function 

global  A  B  G  Z  FUNCNAME  %Matrices  to  be  used 

global  SO  N  BETA  N2  L  CONV  T1  T2%global  variables 

S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-TMTTTT  T  7  AT  TONTS-  S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S'S-S-S-S-S-S-S-S'S-S'S- 

0000000000000000000000-LiN±J.±J_l±ZjiAJ.±WiN  000000000000000000000000000000000000 

SO=0;  %initializing  the  standard  error 

%ll=0;%lf  the  function  is  using  log-likelihood  11==0  least  %squares 
%11==1 

T1  =2 0 ; %input ( ' Enter  the  number  of  iterations?  '); 

N2  =  1 ; %input (' Enter  the  number  of  independent  variables?  '); 

L  =  1 . 0 ; %input (' Enter  the  starting  Lambda?  '); 

CONV  =  0 . 001 ; %input (' Enter  the  convergence  criterion?  '); 

T2  =  0;%actual  number  of  iterations 

s=GET  DATA;%gets  the  data  array  s 
w  =  new  array (s,N2);  %the  data  array 
ydat=w(:,l);  %ydat  gets  Y-data 

for  i=l:N2 

xdat ( : , i ) =w ( : , i+1 ) ;  %xdat  gets  x  data 
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end 

N=length (ydat) ; 

b  num=GET  FUNCTION; %gets  the  function  name  from  the  user 
ll=binary (ydat) ; %binary  determines  if  11=0  or  1 

A=zeros (b  num,b  num) ; %VARIANCE-COVARIANCE  matrix  init.  to  zeros 
G=zeros (b  num, 1 ) ; %gradient  vector  initialized 
modi (b  num, ydat, xdat, 11) ; %begin  parameter  search  modl-modS 
%to  use  these 

%u=mod8 (b_num, ydat, xdat, N-b_num, 11) ; %get  confidence  region  and  %return 
conf.  region  data 

diary  off;  %turns  off  the  diary  file 
return 

%%%%%%%%%%%%%%ENDMARQUARDT%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%BEGIN  MOD  1  %%%%%%%%%%%%%% 
function  modi (b  num, ydat, xdat, 11) 
global  BETA  G  SO  A  N  T2 

global  Z  %Z  is  used  in  mod  6  to  avoid  recomputation  of  p  for  the  SSE  case,  in 
the  LL  case  it  is  not  used  since  we  need  to  recompute  new  partials,  since 
creation  of  an  array  of  dummy  x-values  p  needs  to  be  computed 
p=zeros(N,b  num);  %partial  derivative  in  module  2  is  size=  [N,b_num] 
el=E (ydat, xdat, 11) ;  %function  call  to  get  error  vector 
if  ll==0%ll=test  statement  to  see  if  Log-like  or  least  squares 
SO=sum ( -el ) ; %f or  log-likelihood 
else 

SO=sum (el . ^2 ) ; %f or  least  squares 

end 

T2=T2+1; 
for  i=l : b  num 

BETA ( i ) =BETA ( i )* 1 . 001 ; %BETA  gets  changed  temporarily 

e2 (:,  i) =E (ydat,  xdat,  11)  ;  %e2  is  the  new  error  term  to  be  compared  %to 
el 

BETA ( i ) =BETA ( i )/ 1 . 001 ; %Beta  gets  changed  back 

s ( : , i) = (el-e2 ( : , i) ) ;  %s=temp  array  to  hold  subtracted  values  for  %each 
X  Y  pair,  to  dec.  #  calc 
end 

for  i=l : b  num 

p (:, i ) =s (:, i )/ (BETA ( i )* 0 . 001 ); %creating  the  partial  derivative  for  %each  y 
and  beta 
end 

Z=p;%needed  if  SSE  in  mod  6 
for  i=l : b  num 
if  11==0 

G ( i ) =sum ( -p ( : , i ) ) ; %gradient  vector  for  log-likelihood 
else 

G ( i ) =sum ( (p ( : , i ) . *el ) ) ; %gradient  vector  determining  the  next 

%array  multiplication  is  used  and  not  matrix  multiplication 

end 

end 

for  i=l : b  num 

for  j=l :b  num 

A ( i , j ) =sum (p(:,i) .*p(:,j)); % variance -covariance  vector 
end 

end 

mod2 (b_num, ydat, xdat, 11) ; 
return 

%%%%%%%%%%%%%%END  MOD  1  %%%%%%%%%%%%%% 
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%%%%%%%%%%%%%%BEGINMOD2%%%%%%%%%%%%%% 
function  mod2 (b  num, ydat, xdat, 11) 
global  L  BETA  G~A 
for  i=l ; b  num 

for  j  =1 : b  num 

q(i,j)=A(i,j)/( sqrt ( ( A ( i , i ) *A ( j , j ) ) ) ) ; 

end 

G  (i)  =G  (i)  /  sqrt  (A  (i,  i)  )  ; 

end 

mods (b_num, ydat, xdat, q, 11) ; 
return 

%%%%%%%%%%%%%%END  MOD2%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%BEGIN  MOD3%%%%%%%%%%%%%% 

function  mod3_3 (b_num, ydat, xdat, q, 11) 

global  L  BETA  G  A 

p=zeros ( 1 , b_num) ; 

for  i=l : b  num 

q ( i , i ) =q ( i , i ) * ( 1+L) ;  %adds  the  gradient  vector  to  q 

end 

c=inv(q);  %the  inverse  of  q 

for  i=l : b  num 

for  j  =1 : b  num 

p ( i ) =p ( i ) +c ( i , j ) *G ( j ) ; %estimate  new  partials 

end 

p  (i)  =p  (i)  /  sqrt  (A  (i,  i)  )  ; 

end 

modi (b_num, ydat, xdat, q, p, 11) ; 
return 

%%%%%%%%%%%%%%END  MOD3%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%BEGINMOD4%%%%%%%%%%%%%% 

function  modi (b  num, ydat, xdat, q, p, 11) %q  only  needed  for  recursive  call  %to 
mods 

global  L  BETA  G  SO  A  CONV  T1  T2 
if  T1<0 

disp ( [ ' print  results ' ] ) ; 
mods (b_num, ydat, xdat, 11) ; 
elseif  T1<=T2 %you ' re  out  of  iterations 
tt2=num2str (T2 ) ; 

disp ([' Iteration  number:  '  ,  tt2]); 

T1=T1-1; 

L=0; 

elseif  abs (p . /BETA) <  CONV  %if  any  is  larger  than  conv  do  another  %interation 
disp ( [ ' CONVERGENCE ' ] ) ; 

Tl=-1; 

L=0; 

modi (b_num, ydat, xdat, 11) ; 

else 

BETA=BETA+p; %modif y  parameters  by  adding  their  partial  derivatives 
e=E (ydat, xdat, 11) ; %E  should  return  the  error  for  each  X  Y  pair, 
if  11==0 

sl=sum (-e) ; 
s  s 1 =num2  s  t  r ( - s 1 ) ; 

disp ([' Log-likelihood  for  next  B=  '  ,  ssl]); 
else 

sl=sum (e . ^2 ) ; %summing  up  the  squared  errors 
ssl=num2str (si) ;  %changed  to  positive  from  negative 
disp(['SSE  for  next  B=  '  ,  ssl]); 


end 

if  sl>SO 
L=L*10; 

BETA=BETA-p; %the  estimate  was  worse  and  L  is  incremented  by  then, 

%BETA  is  modified  back  to  its 

%original  value  and  module  3  is  called  again,  this  is  not  counted  %as 
an  iteration 

mod3 (b_num, ydat, xdat, q, 11) ; %Recursive  call  to  mod3 
else 

L=L/10; 

modi (b  num, ydat, xdat, 11) ; %call  modi  again  and  start  over  from  %scratch, 
i.e  new  iteration 

end 

end 

return 

%%%%%%%%%%%%%%END  MOD4%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%BEGINMOD5%%%%%%%%%%%%%% 
function  modS  5 (b  num, ydat, xdat, 11) %not  needed  tl  t2  1 
global  BETA  SO  N  A 

df  =  N-b  num; %the  degrees  of  freedom 
if  11==0~ 

V=l;%for  the  log-like  case 
so=num2str (-SO) ; 

disp([ 'Final  log-likelihood=  ',so]); 
else 

V=SO . /df ; %summing  up  the  squared  errors,  i.e.  variance 

Vl=sqrt (V) ; %stdev 

v=num2str (V) ; 

vl=num2str (VI ) ; 

so=num2str (SO) ; 

disp ( [ ' Variance=  ',v]); 

disp ( [ ' Std . dev=  ',vl]); 

disp([ 'Final  SSE=  ',so]); 

end 

C=inv (A) ; 

A=V*C;  %Var-covar  matrix 

for  i=l : b  num 

D ( i ) =sqrt (A ( i , i ) ) ; %D=std .  error  of  parameter 

end 

D2=D./BETA;  %coefficient  of  variation 
b=num2str (BETA) ; 
d=num2str (D) ; 
d2=num2str (D2 )  ; 

disp ( [ ' Parameters=  ' , b] ) ; 

disp(['Std.  Error  of  Parameters=  ',d]); 

disp(['Coeff  of  var=  ',d2]) 

disp ( [ ' Var-Covar  matrix=  ' ] ) ; 
disp (A) 
return 

%%%%%%%%%%%%%%END  MOD5%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%BEGINGET_FUNCTION%%%%%%%%%%%%%% 
function  y=GET_FUNCTION 

global  FUNCNAME  npars  BETA%FUNCNAME  is  global  and  used  by  other  %functions  to 
call 

%the  chosen  function,  contains  the  anme  of  the  function  and  its  %destination 
%  select  model 

[funcfilename,  pathname]  =  ... 
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uigetfile ( 'AB_* .m' Select  model  function  file'); 

FUNCNAME  =  strtok ( f uncf ilename , 

%THE  BELOW  COMMAND  IS  HOW  TO  EVALUATE  THE  FUNCTION  I.E.  CALL  FUNCNAME 
%string  =  [ ' u= '  funcname  ' (ydat, W, BETA) ; ' ] ; %creates  a  string  to  be  used  %in 
eval 

%eval ( string) ;  %the  driver  to  call  the  function 

Z  =  [le-10  le-10  le-10  le-10  le-10  le-10  le-10  le-10  le-10  le-10  ] ; 

string  =  [ ' dummy= '  FUNCNAME  '  ( 0 , 0 , Z ) ;  ' ] ; 

eval ( string) ;  %  dummy  call  to  funcc  to  get  npars 

clear  Z; 

for  i  =  l:npars%get  intial  betas  from  user 
num  param  =  int2str(i); 

ask  value  =  ['Enter  initial  value  for  beta ( '  num  param  '):  ' ] ; 

BETA(i)  =  input (ask  value); 

end 

y=npars ; 
return 

%%%%%%%%%%%%%%ENDGET_FUNCTION  %%%%%%%%%%%%%% 
%%%%%%%%%%%%%%BEGINGET_DATA%%%%%%%%%%%%%% 
function  y=GET_DATA 

disp(' Please  provide  the  name  of  for  the  array  of  outcome  data.'); 

[filename, pathname] =uigetfile ('*. txt ',' Outcome  Data  File  Name ' ) ; %gives  %the 
name  of  file  and  path 

disp(['Data  file  name:  '  filename] ); %displays  the  file  name  chosen 
fullfilename= [pathname  filename] ; 

%s=dlmread (fullfilename, ' , ' ) ; %reads  the  data  array 

in  =  fopen (fullfilename, ' rt '); %fopen (filename, permission) ,  rt=rad  and  %write 
a  txt  file 

s  =  fscanf(in, ' %f ' ) ; %reads  the  entire  file  into  an  array,  file  needs  to  %have 

the  y  and  x  variables  column  wise 

%i.e.Y,  X(l),  X(2),  X(3),  ...X(N2)  etc 

y=s;%returns  s,  i.e.  the  whole  data  array 

fclose(in);  %closes  data  file 

return 

%%%%%%%%%%%%%%END  GET_DATA%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%BEGIN  E  %%%%%%%%%%%%%% 
function  y=E (ydat, xdat, 11) 
global  FUNCNAME 
global  BETA 

str  =  [ ' v= '  FUNCNAME  ' ( ydat , xdat , BETA) ;']; %creates  a  string  to  be  used  %in 
eval 

eval (str);  %the  driver  to  call  the  function 
for  i=l : length (ydat) %to  be  used  when  v  cannot  be  0 
if  v(i)==0 

v(i)=0. 000001; 
elseif  v(i)==l 

v(i)==0. 999999; 

end 

end 

if  11==0  %if  LL=compute  log-likelihood  for  each  observation 
m=ydat . *log (v) + ( 1 . 0-ydat) . *log ( 1 . 0-v) ; 
else  %if  SSE  return  the  error  estimate 
m=ydat-v; 

end 

y=m; %returns  LL  or  SSE  also  called  F  in  all  NMRI  programs 

%%%%%%%%%%%%%%END  E%%%%%%%%%%%%%% 
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%%%%%%%%%%%%%%BEGINNEW_ARRAY%%%%%%%%%%%%%% 

function  y=  new  array (W,N2) 

v=0;  %coounter  in  the  for  loop 

col=N2+l;  %col  gets  total  number  of  columns 

row=length (W) /col ; %this  is  the  number  of  data  points  in  each  column,  %i.e. 
rows 

for  i=l ; row 

for  k=l : col 
v=v+ 1 ; 

temp ( i , k) =W (v)  ; 

end 

end 

y=temp; % returns  the  reshaped  matrix  with  Y  in  the  first  column  and  then  %the 

X '  es 

return 

%%%%%%%%%%%%%%END  NEW_ARRAY%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%BEGINMOD8%%%%%%%%%%%%%% 
function  y=mod8 (b_num, ydat,  xdat,  df ,  11) 

global  A  N  Z  BETA%Z  is  a  matrix  containing  the  partial  derivative  p()  from 
modi 

T=l/ (df) ; 

T=1.96+T* (2.3724+T*(2.8227+T*(2.5561+T*1.5897))); 
if  df<=l.l 
T=12.706 

end 

xmin  =  min (xdat);  %gets  minimum  x 

xmax  =  max (xdat);  %gets  max  x 

xrange  =  xmax  -  xmin;  %gets  the  range  of  values 


if  ll==0%if  LL=set  all  outcome  to  1  and  reestimate  partials  based  on  this 
newx  =  linspace (xmin, xmax, 1 00+1 ) ' ; %creates  a  data  array  of  values  at  each 

who 

Vl=zeros (length (newx) , 1) ; %initializing 
YHAT=zeros (length (newx) , 1) ; %initializing 

yones=ones ( length (newx) , 1 ); %creates  an  array  with  all  ones  for  use  in 
making  conf.  region 

p=zeros ( length (newx) , b  num) ;  %partial  derivative  in  module  2  is  by  size= 
[N,  b_num] 

el=E (yones,  newx,  11) ;  %function  call  to  get  error  vector  or  LL  also  F  in 
the  NMRI  programs 
for  i=l : b  num 

BETA ( i ) =BETA ( i )* 1 . 001 ; %BETA  gets  changed  temporarily 

e2 (:, i) =E (yones, newx, 11) ;  %e2  is  the  new  error  term  to  be  compared  to 
El 

BETA ( i ) =BETA ( i )/ 1 . 001 ; %Beta  gets  changed  back 

s ( : , i) = (el-e2 ( : , i) ) ;  %s=temp  array  to  hold  subtracted  values  for  each 
X  Y  pair,  to  dec.  #  calc 
end 

for  i=l : b  num 

p  (:, i ) =s (:, i )/ (BETA ( i )* 0 . 001 ); %creating  the  partial  derivative  for  each 
y  and  beta 
end 

for  q=l : length (newx) 
for  i=l : b  num 

for  j=l :b  num 

VI (q) =V1 (q)+(p(q,i) .*p(q,j) .*A(i,j)); % variance -covariance  vector 
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end 


end 

end 

SEYHAT=sqrt (abs (VI ) ) ;  %the  seyhat  of  the  estimate 


YHAT=exp (el ) ;  %for  LL  estimation 

err  maxmin=T*SEYHAT; %for  LL  estimation 

YMIN=max (0, exp (el-err  maxmin) ) ; %The  lower  end  of  the  95%CL  region 
YMAX=min ( 1 , exp (el+err  maxmin) ); %The  upper  end  of  the  95%CL  region 


else  %if  the  estimate  is  continuous  variables  and  uses  SSE  estimation 
Vl=zeros (length (xdat) , 1 ) ; %initializing 
YHAT=zeros (length (xdat) , 1 ) ; %initializing 

el=E (ydat, xdat,  11) ;  %function  call  to  get  error  vector  or  LL  also  F  in  the 
NMRI  programs 

for  q=l : length (xdat) 
for  i=l : b  num 

for  j  =1 : b  num 

VI (q) =V1 (q)+(Z(q,i) .*Z(q,j) .*A(i,j)); % variance -covariance  vector, 
uses  Z  vector  from  modi 
end 

end 

end 

SEYHAT=sqrt (abs (VI ) ) ;  %the  seyhat  of  the  estimate 
YHAT=ydat-el ;  %as  defined  in  Bailey  and  Homer 
YMAX=YHAT+T* SEYHAT;  %for  SSE  estimation 
YMIN=YHAT-T* SEYHAT; 
newx=xdat; 

end 


IF  ELSE  STATEMENT++++++++++++++++++++++++++ 


plot (newx, YHAT, ' b : p ' , newx, YMIN, ' c- ' , newx, YMAX, ' c- ' ) 
disp(['X  YHAT,  SEYHAT,  YMAX,  YMIN']) 

u=[newx,  YHAT,  SEYHAT,  YMAX,  YMIN] ; 
u=num2str (u) ; 
disp (u) 

y=u;%return  data  to  marquardt 
return 

%%%%%%%%%%%%%%END  MOD8%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%BEGINBINARY%%%%%%%%%%%%%% 
function  y=binary (bin) 
for  i=l : length (bin) 

if  (bin(i)==l)  |  (bin ( i ) ==0 ) %if  outcome  1  or  zero,  i.e.  binary 

y=0; 
else 

y=l;  %if  each  data  not  binary  it  use  SSE  and  set  11=1 

return  %if  not  1  and  0  return  immediatley  no  sense  to  continue 

end 

end 

return 

%%%%%%%%%%%%%%END  BINARY%%%%%%%%%%%%%% 
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Variables  in  Use: 


A(N,N)=holds  information  for  variance-covariance  matrix 
G(N)=gradient  vector 

FUNCNAME  =  string  for  name  of  function 
SO  =  sum-squared  errors 
N  ^number  of  observations 
N1  =  number  of  parameters 
N2  =  number  of  independent  variables  (x) 
b  num  =  number  of  parameters 
BETA  =  array  of  parameters 
L  =  Lambda 

CONV  ^convergence  criterion 
T 1  =  maximum  number  of  iterations  before  stopping 
T2  =  iteration  counter 

V  =  variance  of  observation 

V I  =  standard  deviation 

11  =  holds  the  type  of  dependent  variables,  binary  or  continuous 

p(N,  b  num)  =  partial  derivatives 

xdat  =  independent  variable  array 

ydat  =  dependent  variable  array 

e  =  error 

e  1  =  error 

e2  =  error 

Z  (N,  b  num)  =  temporary  global  storage  of  partial  derivative 
q  =  intermediate  storage 
q,  i,  j  =  counter  variables 
df  =  degrees  of  freedom 
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APPENDIX  B 


To  run  the  data  in  Table  1  using  the  non-linear  least  squares  method,  Eq.  1  will  be  used  to 
fit  the  data.  First  the  data  are  saved  as  a  text  file  and  a  Matlab  m.file  is  made  that  defines  the 
equation.  In  this  case,  the  function  file  is  called  AB  Homerl: 

function  y=AB  Homer 1 (ydat, xdat, BETA) 
global  npars 

if  BETA(l) ==le-10  &  BETA (2 ) ==le-l 0 
npars  =  2 ; 
return 
else 

y  =  BETA(l)  *  (xdat  (  :  ,  1)  . (BETA(2)  *xdat  (  :  ,  2)  )  )  ; 

end 

return 

The  function  can  be  called  any  name,  beginning  with  “AB_”.  The  “AB_”  can  be  changed  by 
modification  of  the  “Get_Function”  m.file.  Next,  “N2”  in  the  “marquardt”  routine  is  set  to  2  for 
two  independent  variables,  convergence  (CONV)  to  0.001,  Fambda  (F)  to  1,  and  the  program 
run  with  a  maximum  of  20  iterations  (Tl).  Starting  values  for  the  parameters  are  2  and  2  for  pi 
and  P2,  respectively.  The  output  of  the  result  is  shown  below  and  is  equivalent  to  the  result 
presented  by  Homer  and  Bailey  (1). 

Enter  initial  value  for  beta(l):  2 
Enter  initial  value  for  beta(2):  2 
SSE  for  next  B=  4.9452 
SSE  for  next  B=  4.7942 
SSE  for  next  B=  4.7503 
SSE  for  next  B=  4.75 
CONVERGENCE 
print  results 
Variance=  1.1875 
Std.dev=  1.0897 
Final  SSE=  4.75 

Parameters^  2.2499  1.8301 

Std.  Error  of  Parameters^  0.54486  0.4006 

Coeffofvar=  0.24217  0.21889 

Var-Covar  matrix^ 

0.296875  -0.190239 
-0.190239  0.160477 
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APPENDIX  C 


The  maximum  likelihood  estimation  procedure  is  used  to  fit  the  data  in  Table  2.  The 
outcome  of  the  data  is  binary,  and  in  this  case  a  response=l,  and  no-response=  0.  The  data  is  fit 
to  the  dose  response  functions  Eq.  2  and  Eq.  3. 

Example  1 

Eor  Eq.  2,  the  parameter  (Pl)  is  the  dose  at  which  50%  of  the  outcome  is  1.  The  function 
file  for  this  equation  is  called  AB_Homer2  ; 

function  y=AB  Homer2 (ydat, xdat, BETA) 
global  npars 

if  BETA(l) ==le-10  &  BETA (2 ) ==le-10 
npars  =  1; 
return 
else 

y  =  xdat ./ (BETA ( 1 ) +xdat)  ; 

end 


Eor  this  estimation,  the  “N2”  in  “marquardt”  needs  to  be  set  to  1,  while  all  other  variables  remain 
as  in  the  example  above.  The  output  of  the  results  for  a  starting  value  of  the  parameter  of  0.45  is 
shown  below: 

Enter  initial  value  for  beta(l):  0.45 
Log-likelihood  for  next  B=  -6.8292 
Log-likelihood  for  next  B=  -6.827 
Log-likelihood  for  next  B=  -6.827 
Log-likelihood  for  next  B=  -6.827 
CONVERGENCE 
print  results 

Final  log-likelihood=  -6.827 
Parameters^  0.49645 

Std.  Error  of  Parameters^  0.33669 
Coeffofvar=  0.67819 

Var-Covar  matrix= 

0.113360 

Example  2 

Equation  3,  also  known  as  the  Hill  Eunction,  is  commonly  used  to  describe  biological 
phenomena  such  as  the  02-dissociation  curve  (2),  and  has  been  used  to  describe  the  probability 
in  decompression  sickness  (7,8).  Eor  our  purpose,  we  use  it  to  describe  a  dose-response 
relationship  with  only  two  outcomes.  Again,  the  first  parameter  (Pi)  is  the  dose  at  which  50%  of 
the  outcome  is  1,  and  the  second  parameter  (P2)  is  the  slope  of  the  sigmoidal  dose  response 
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curve.  Compared  to  AB_Homer2,  the  funetion  file  for  this  equation  only  requires  the  following 
ehanges: 

npars  =  1; 

and 

y  =  xdat . / (BETA ( 1 ) +xdat) ; 

The  “marquardt”  remains  the  same  as  for  the  first  example  using  maximum  likelihood  above. 
The  printout  of  the  result  is  shown  below  using  the  following  starting  parameters  Pi  =0.8613  and 

Pi2=3.5338,  and  with  Lambda=  1.0,  eonv=0.001,  andTl=20: 

Enter  initial  value  for  beta(l):  0.8613 
Enter  initial  value  for  beta(2):  3.5338 
Log-likelihood  for  next  -5.5912 

Log-likelihood  for  next  -5.5904 

Log-likelihood  for  next  B=  -5.5902 
Log-likelihood  for  next  B=  -5.5906 
Log-likelihood  for  next  B=  -5.5905 
Log-likelihood  for  next  B=  -5.5902 
Log-likelihood  for  next  B=  -5.5902 
Log-likelihood  for  next  B=  -5.59 
Log-likelihood  for  next  B=  -5.5904 
Log-likelihood  for  next  B=  -5.5901 
Log-likelihood  for  next  B=  -5.59 
Log-likelihood  for  next  B=  -5.59 
Log-likelihood  for  next  B=  -5.59 
Log-likelihood  for  next  B=  -5.59 
CONVERGENCE 
print  results 

Final  log-hkehhood=  -5.59 
Parameters^  0.86235  3.4309 

Std.  Error  of  Parameters^  0.34462  2.005 

Coeffofvar=  0.39962  0.58439 

Var-Covar  matrix^ 

0.11875  0.55232 
0.55232  4.02002 
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